'''
Created on May 21, 2011

@author: oabalbin
'''

def read_gene_list(gfile):
    '''
    '''
    gfile = open(gfile)
    glist=[]
    for l in gfile:
        glist.append(l.strip('\n'))
    return glist

def read_cons_table(ifile,number_Alg, glist, ofile):
    '''
    '''
    
    SNPAlg=number_Alg
    Major=1
    SNP_ANNOT=3 # Actionable, HapMap, dbSNP
    NFIELDS=9
    VARS=1
    SAM=0
    GATK=1
    IMAT=SNPAlg+Major+SNP_ANNOT
    # SAM
    M0_QUAl=IMAT + NFIELDS*SAM + 1 #14 sam
    # GATK
    M1_VQSLOD=IMAT + NFIELDS*GATK + 8
    M1_QUAl=IMAT + NFIELDS*GATK + 1
    M1_DP=IMAT + NFIELDS*GATK + 2
    M1_SB=IMAT + NFIELDS*GATK + 7
    #
    ###
    mut_type=IMAT + 2*NFIELDS*GATK+1
    gene_name =IMAT + 2*NFIELDS*GATK+2
    MajCol=3
    VarSCol=2
    SamCol=1
    GatkCol=2
    Majority=2
    Mone=1
    Mtwo=2
    gpoints=1000
    ipoints=500
    TRUE_SNPs= IMAT #IMAT=dbSNP, IMAT-1=Hapmap
    
    thismuts=5.0
    dbsnpth=1.0
    DP_Threshold=10.0
    SB_Threshold= -0.01
    
    ifile=open(ifile)
    header = ifile.next().strip('\n').split('\t')
    ofile = open(ofile,'w')
    header.append('Actionable2')
    print header
    ofile.write(",".join(header).replace(',','\t')+'\n')
    for l in ifile:
        fields = l.strip('\n').split('\t')
        # filter non-synonymous and not in dbSNP snvs
        # Filter based on being called for both or algorithms or gatk
        
        #print header[mut_type], mut_type, IMAT, NFIELDS*GATK
        #print header[MajCol],header[GatkCol] 
        if ((float(fields[mut_type]) <= thismuts) and (float(fields[TRUE_SNPs]) < dbsnpth)) and \
           (float(fields[MajCol]) == Majority or float(fields[GatkCol]) == Mone):       
            #total depth of coverage and strand bias
            if (float(fields[M1_DP]) >= DP_Threshold and float(fields[M1_SB]) <= SB_Threshold ):
                thg = fields[gene_name].split('(')[0].strip('\'')
                if thg in glist:
                    fields.append('1')
                    print fields[gene_name]
                else:
                    fields.append('0')
                
                ofile.write(",".join(fields).replace(',','\t')+'\n')
            
gfile='/exds/users/oabalbin/projects/exome/actionable_gene_list.txt'        
ifile='/exds/projects/exome/FMPN_51/analysis/cons_calls/concesus_table_FMPN51_notBenign.txt'
ofile='/exds/projects/exome/FMPN_51/analysis/cons_calls/concesus_table_FMPN51_notBenign_report.txt'

number_Alg=2
genelist = read_gene_list(gfile)
read_cons_table(ifile,number_Alg,genelist, ofile)